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We systematically calculate thermopower of biased and unbiased multilayer graphene systems. 
The effect of screening to a bias field perpendicular to the graphene planes is taken into account self- 
consistently under the Hartree approximation. The model including nearest neighbor hopping and 
the more general Slonczewski-Weiss-McClure (SWMcC) model are both considered for a comparison. 
The effect of impurity scattering is studied for monolayer and unbiased bilayer graphene and is 
treated in terms of the self-consistent Born approximation. For a monolayer graphene, only when 
the effect of impurity scattering is taken into account, could all the qualitative aspects of the 
experimental results be correctly reproduced. Besides bilayer graphene, only trilayer graphene 
opens a small gap and shows a slight enhancement of thermopower under an external bias. The 
biased bilayer graphene shows the largest thermopower among all the systems studied. 

PACS numbers: 72.80.Vp, 72.10.-d, 73.50.Lw 
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I. INTRODUCTION 

Recently, thermopower of graphene systems have at- 
tracted much attention. Il3jj In experiments on mono- 
layer graphene (MLG), while the high carrier density re- 
gion is well accounted for in terms of the Mott's rela- 
tionship, the low carrier density region shows deviation 
away from this behavior. 0-0] This is explained in terms 
of either an electron-hole puddle model [1] or the coher- 
ence effect between the conduction band and the valence 
band mediated by impurity scattering 0- Besides the 
MLG, a bilayer graphene (BLG) system is also very in- 
teresting. The opening of a gap by an external electric 
field applied perpendicular to the layer planes [l3 - il9| in- 
troduces a new degree of tunability to the system and 
is recently shown by the authors to enhance greatly the 
thermopower . 

Despite the above progresses, several questions remain 
to be answered. First, the dependence of thermopower 
on carrier density as a function of temperature for a MLG 
has two features which are not fully explained in previous 
works. The first one is the shift of thermopower peak 
with respect to the charge neutral point as temperature 
decreases. In addition, thermopower shows a monotonic 
dependence on temperature for carrier densities smaller 
than that of the thermopower peak.0, 0, !| The second 
question is the systematic dependence of thermopower 
on the number of layers. Further more we are interested 
in finding out whether or not an applied biased potential 
will enhance thermopower for the multilayer graphene as 
the same as the bilayer system. Below we shall address 
these questions. 

Full microscopic calculations of thermopower of MLG 
in the presence of charged impurity scattering are per- 
formed in terms of the self-consistent Born approxima- 
tion. Better agreement with experiments @tJ] as com- 
pared with previously reported results are achieved. 0, 
ID, l2( We also study thermopower of impure unbiased bi- 
layer graphene. Results which are qualitatively similar to 
those in monolayer graphene are obtained, but for much 



higher impurity concentrations. 

For a multilayer graphene with layer number from 
two to six, we consider both the nearest neighbor tight 
binding model and the mo re g eneral Slonczewski-Weiss- 
McClure (SWMcC) model. HIH For a multilayer sys- 
tem the effect of screening is important in determining 
the band structure, as it produces non-uniform charge 
densities across the layers. [13, EH HHIll In this work, 
by using the self-consistent Hartree approximation we 
obtain potential energies of different layers for unbiased 
systems with three or more layers and biased systems 
with a perpendicular electric field for two or more layers. 
For unbiased graphene multilayers, peak values of the 
thermopower are close to each other irrespective of the 
model used. For a biased multilayer graphene, only ther- 
mopower of the bilayer shows significant enhancement, 
thermopower of the trilayers also increases slightly, but 
thermopower of systems with more layers decreases un- 
der a bias. The huge enhancement of thermopower in a 
bilayer graphene is a direct result of gap opening under 
a perpendicular electric bias.[ll| In a trilayer graphene, 
a much smaller gap opens and hence the thermopower 
enhancement is very small. While for layer number 
larger than three, no gap opens under a bias. Thus no 
thermopower enhancement is expected for a multilayer 
graphene with layer number larger than three. 



II. MODEL AND METHOD 

A. The SWMcC model 

We consider multilayer graphene systems stacked in 
the standard AB (Bernal) stacking between consecutive 
layers. |29rf33 | In order to correctly reproduce the band 
structure, a tight binding type of model up to the fifth 
nearest neighbor hopping is employed. [20l - [22l l29l. l3ll HH 
For an iV-layer graphene, the Hamiltonian is transformed 
in the xy-jAane to wave vector space, whereas the direc- 
tion perpendicular to the plane is kept in the real space. 
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Along the z-direction, the lattice sites which are verti- 
cally aligned to each other are labeled as (Ai, B 2 , A 3 , 
...), with Ai (Bi) denoting the A (B) sublattice atoms 
on the i-th lay er. 13 Cll The single spin Hamiltonian is thus 
written asfH |29U3ll 132^ 



N 
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n , k 



+72(^2 n-l,k^2n-l,k + a 2n,k a 2",k)] 

+ X][ 7l ( a 2n-l,k & 2n,k + aL+l,k & 2n,k + Ha 



, 75 / f 

+ yl a 2n-l,k a 2n+l:k 



6 2n,kW2,k + #-C.) 



12 + t 

+ Y ( & 2n-l,k 6 2n+l,k + <4 n k a2„ +2 ,k + H.C.)} 

+ v 3 5Z^(k)&L-l,k a 2r 1 ,k + ^(k)4n+l,k a 2 n ,k + H.C.} 



-Vi 



5Z[^*( k )4n-l,k a 2n,k + 0*(k)& 2 „-l : k & 2«,k 



+0* (k)4 n+1 k o 2n ,k + 4>* (k)4„ + i.k 6 2«.k + H.c] . (1) 



in which 



ffo(k) = 





0*(k) / 



0(k)^ (0 71 







(4) 



In the absence of any external field or gate voltage, 
it is known that the above Hamiltonian could be de- 
composed into subsystems of bilayer and monolayer 
graphene. [3(1 IH, HI, H3 For an odd-layered multilayer 
graphene, there are one MLG subsystem and (N — l)/2 
different BLG subsystems. While for an even-layered 
multilayer, there are N/2 different bilayers. Label the 
different BLG subsystems with an index to, and take the 
basis vector as Vi(k)=(4a k , 6^ lk , a^ 2k , b^), the 
Hamiltonian matrix of the m-th BLG subsystem is writ- 
ten as 



( <f>(k) 7 i„ 

»*(k) 

cj)(k) 

\ 7lm 0*(k) 



(5) 



where 7i m =27i sin 
N — 3, 



takes the value of N — 1 , 



2(iV+l) ' 

, 2 (for odd N) or 1 (for even N) 



C. Screening effects 



In the above expression, 0(k)=7o exp(ik • di) arises 

1=1 

from motion within separate layers. ^3=73/70, 
^4=74/70- 7o, 7i, 72, 73, 74, 75 and A are standard 
SWMcC parameters, and are taken as 3.12 eV, 0.377 
eV, -0.02 eV, 0.29 eV, 0.12 eV, 0.0125 eV and 0.004 
eV, respectively. [i3,|2l|2i,|2ll lag Except the first term, 
summation over n runs from 1 to N/2 or N/2 ± 1 etc., so 

that the basis vector ip" 1 (k) = (a[ k , b\ k ; a^, b\ k ;...; a Nk> 

^Ivk) nas 2iV components. The Hamiltonian is simply 
given as 



(2) 



where matrix form of (k)is not shown explicitly here 
but can be found in Refs. |22f and 1321. 



B. Nearest neighbor tight binding model 

Apart from the full SWMcC model, a simplified model 
up to nearest neighbor intralayer and interlayer hopping 
is widely used to study many problems. In this case, the 
Hamiltonian matrix is written as 129, •'!(). liA Kill 



ff(k) 



/ifo(k) V 
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V H {k) V 
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In our former work, a bias voltage is shown to greatly 
enhance the thermopower of BLG. 11] It is unclear 
whether this is true also for other multilayer graphene. 
In the presence of a bias voltage or a charged gate, 
the potential energy difference would be induced be- 
tween different layers. To determine the distribution 
of potential energies, self-consistent calculations are re- 
quired to take into account of the Coulomb screening 
arising from char ge redistribution among the different 
layers. [iH I22I I24I |38| For unbiased multilayer graphene 
with layer number larger than two, the charge and po- 
tential energy distribution should also be obtained self- 
consistently. 

There are two schemes to account for the charge redis- 
tribution. One is to consider the screening of a charged 
gate with fixed carrier density by the multilayer graphene 
system 22] . The other approach is to consider the charge 
redistribution problem of the multilayer graphene system 
in the presence of a perpendicular electric field, while the 
carrier density is controlled by another gate voltage ap- 
plied equally to all layers [38j. Here we follow the second 
approach. 

Suppose an external electric field E is applied to the 
free-standing multilayer graphene system along the di- 
rection perpendicular to the graphene planes from layer 
1 to layer N (Eo=0 for unbiased systems). After estab- 
lishing equilibrium, the i-th layer has a total excess elec- 
tron density of nj. According to Gauss's law, the screen- 
ing electric field pointing from the j-th to the 0+1)- 

N j 

th layer is E s ^ j+1) ^^-( m ~ S where 'sc' 

i=i+l 



i' = l 
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means screening, e is the absolute value of the electron 
charge, and e r =3 is the static dielectric constant of the 
multilayer graphene [3] . The total electric field is thus 
E(j j +1 }=Eo + Efjj+iy Introduce the doping Xi=riiQo 

(Qo=^j^a 2 is the effective area per carbon site), which 
means the number of excess electrons per site residing 
on the i-th layer. Define two parameters Vo=eEodo and 
7= 2 e ^ . do— 3.5 A is the vertical distance between adja- 
cent fayers. The potential energy difference between the 
j-th and the (j+l)-th layer is 

N j 

AVj = V j+1 - Vj = V + 7 [ Yl x i ~ E ( 6 ) 

i—j + l i' — 1 

The total difference of on-site potential energy is thus 



JV-l 



v N -v l = Y J 



JV-l 



= (N - 1)V + 7 E ( N ~ Mxn-j+i - xj). (7) 

Since only differences in the potential energies are rele- 
vant, we would take the potential energy of the first and 
N-th layer symmetric with respect to zero by shifting 
the energy reference. Thus we take Vm—{Vn — Vi)/2 and 
V\=-Vn- Potential energy differences Vj+i — Vj between 
consecutive layers are not influenced by this shift of en- 
ergy reference, and are still determined by Eq. ([5]). 

For a certain Vq and a fixed average doping, we cal- 
culate the set of doping {a;,;} self-consistently. Label 
the 2Nx2N matrix which diagonalizes the Hamilto- 
nian matrix as U(k), that is [ft(k)if (k)[/(k)=ff d (k). 
Here, we have added the on-site potential energies 
Hy=J2 Vv^vO c )'0k to the model and have kept the no- 

k 

tation of the Hamiltonian matrix unchanged as H(k). 
The diagonal matrix Hy(k) is k independent and is de- 
fined as diag[t^i, Vi, ... , Vn, Vn]- The a-th column 
of U(k) stores the a-th eigenvector of H(k) correspond- 
ing to an eigenenergy of e ct (k)=[i7 ( i(k)] aQ . Suppose the 
number of wave vectors considered in the Brillouin zone 
is iVk, the self-consistency condition for determining n 
is thus 

Xi + 1 = J£ Ed^-^^l 2 + l^ Q (k)| 2 )/(e«(k) - M ). 

k , a 

(8) 

[i is the chemical potential. f(x)=l/(e /3x + 1) is the 
Fermi distribution function, in which f3 represents the 
inverse temperature l/fc^T with ks denoting the Boltz- 
mann constant. After a set of convergent results for 
{xi,Vi;i — 1,...,N} is obtained, the corresponding on 
site potential energies Vi (i = 1,...,N) are substituted 
into the model to calculate the thermopower of interest 
in this work. 



D. Scattering by charged impurities 

Transport properties of graphene systems sensitively 
depend on the nature and concentration of impurities. In 
particular, charged impurities are shown to dominate the 
various transport properties of MLG [39l - l42T | and also play 
a very important role in the transport of BLG[43]. In the 
following, we would consider the effect of charged impu- 
rity scattering in the MLG and unbiased BLG systems. 
The treatment would be at the level of self-consistent 
Born approximation (SCBA).@, H, Hi-E^ Due to the 
limitation in computation time and the fact that calcu- 
lations performed for clean systems are enough to an- 
swer if a bias could enhance thermopower of a multilayer 
graphene, in this work we would concentrate on clean 
systems for layer number larger than two. 

The formalism is presented in detail in our former work 
on thermopower of gapped BLGTlj. The thermopower 
is represented in terms of the linear response coefficients 
as[T], 0, EH, 



5* = 



L 



12 



eTL u 



(9) 



where T is the absolute temperature and e is the magni- 
tude of electron charge. The linear response coefficients 
are obtained as 



lim Rc£y (w + iO + ). 



(10) 



In the Matsubara notation, the correlation function 
reads Ei 



iT 



dre^iT^ir)-^)), (11) 



The two subindices 'i' and 'j' both run over 1 and 2. ri is 
the total effective area of the system defined as the area 
per layer multiplied by N, the number of layers. d=2 is 
the dimensionality, ji and j2 denote par ticle current and 
heat current operators, respectively. [1 1| The two linear 
response coefficients in terms of which the thermopower 
is obtained are written as 

Ln= T / +D °^[-Mfl]Re{P lj ( e - J 0+,e + zO+) 



2tt l de 



(12) 



where the subindex j is either 1 or 2. The kernels are 
defined as [HI 

P l3 (z, z') - ]T Tr{G k (z)r!(k, z, z')G^{z') ■ jf }. 

(13) 

with ri(k, z, z') as the vertex function corresponding to 
the wave vector k. Here, z=e ± i0 + and z'—e + i0 + . is 
the matrix for the charge current at wave vector k and 
is obtained from the Hamiltonian matrix -ff(k) of the 
multilayer graphene as [ill US HH 



Ji 



k - V k ff(k). 



(14) 
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In -ff(k) the potential energy distribution Hy(k) with 
the on site energies determined self-consistently are in- 
corporated, which affects the energy spectrum but does 
not affect the current operator. 

For an N-layei multilayer graphene, the Green's func- 
tion is represented in terms of a 2Nx2N matrix as 

G k (z) = [Gi(z)- 1 -^(z)]-\ (15) 

The free Green's function is obtained as G£(z) = [(z + 
^)I2N — -^(k)] -1 , where I2N is a 27V-dimensional unit 
matrix and /1 is the chemical potential determined from 
-ff(k) for a certain carrier density and temperature. 
The self energies and the vertex functions are deter- 
mined self-consistently in terms of the following iteration 
functions [7LITTI] 

Sk W = S^ k(k_k,)|2[G ^ (zrl_Ek ' (z)r1 ' (16) 

k' 



r,(k,3,s')=jf 
♦SEN*- 

k' 



(17) 



k')| 2 G k ,(z)r 1 (k',2 1 2')G k K^)- 



rii is the impurity concentration averaged to per layer. 
Wj(q) is the electron-impurity scattering potential. For 
charged impurity scatte ring , Wj(q) is taken as of the 
Thomas-Fermi type[H, 0, Ell 



Vi(q) 



2ne 2 



£r(<? + QTf) 



(18) 



e r is the effective dielectric constant from lattice and sub- 
strate, e r =3 is adopted in this work@, H^. d{ is the 
distance between the impurities and the graphene plane 
and would be set as zero in the present work [J, Q. qxF 
is the Thomas-Fermi wave number and is obtained from 
the long-wavelength-limit static polarizability of the cor- 
responding noninter acting electron system [IlUfil as 



QTF = 2?re xfcr, 
with the static polarizability 

X=^ f dr(T r n(T)nt(0)) c . 



(19) 



(20) 



The particle number operator is defined as 
n(r)=J3 V ; k( T )'^' k ( T )- A factor of '2' comes from 

k 

the two fold degeneracy in spin. The subindex 'c' 
means retaining only connected Feynman diagrams in 
evaluating the expectation value. [Til] 

The self-consistent Born approximation (SCBA) en- 
ters when making averages over impurity configurations 
as 52] 



(Pi(q)ft(-q')> = N A 



(21) 



where iVj = UiQ is number of impurities in the system 
under consideration. 

For clean systems, Gk(z) reduces to the free Green's 
function and the vertex function r*i ( k, z, z') reduces to 
the bare charge current matrix i k .[ll| 



III. RESULT AND DISCUSSION 
A. Monolayer graphene 

The experimental results 0-Q for thermopower of a 
MLG show two interesting features. First, for a certain 
temperature, thermopower follows 1 / \/\x\ (x is the aver- 
age number of excess electrons per site) for high carrier 
densities but then the magnitude decreases and changes 
sign close to the charge neutrality point. Theoretical 
calculations in terms of both the Boltzmann transport 
equationjlj and the microscopic Kubo's formulaQ have 
successfully reproduced the above features. The devia- 
tion of thermopower from the 1 / y/\x\ behavior close to 
the charge neutrality point is ascribed to electron-puddle 
formation [lj or coherence between the conduction band 
and valence bandQ, which both imply the coexistence of 
electron-like and hole-like characters of the carriers. The 
second feature is that, as temperature decreases the ther- 
mopower decreases monotonously for all carrier densities 
and the peak position also shifts toward lower carrier 
densities. [H, Q The above behaviors are very similar to 
the results presented in Figs. 1(c) and 1(d). 

Results in Fig. 1 are obtained in terms of the near- 
est neighbor hopping model with hopping integral as 2.7 
eV. As shown in Fig. 1(a) and 1(b), calculations on 
clean systems are unable to reproduce the second feature 
mentioned above. [9j Previous works on transport proper- 
ties of a MLG had confirmed the importance of charged 
impurities. [39M42I ] Here, we treat the impurity scattering 
in terms of SCBA.0,E3 We would neglect the shift of 



chemical potential by the impurity potential. [7J, lll| The 
results for a series of temperatures with impurity concen- 
tration Tii — 10 12 cm~ 2 , which corresponds to a moder- 
ately disordered sample |40j, are shown in Fig. 1(c). For 
the electron-hole symmetric band as considered here, the 
relationship S{— x) = —S(x) generally holds. The full 
curve for 300 K in Fig. 1(c) is explicitly calculated, the 
above relationship is seen to hold very well. For all other 
curves, only the hole doping part with x < is calculated 
explicitly and the part with x > is obtained in terms 
of S(x) — —S(—x). It is clear that, Fig. 1(c) agrees 
qualitatively very well with experiment by Zuev et al. 0] 
As shown in Fig. 1(d), increasing m to 5 x 10 12 cm~ 2 , 
which corresponds to dirty graphene samples, the peak 
positions are shifted to higher carrier densities and the 
maximum values of thermopower are further suppressed 
for all temperatures, which are qualitatively similar to 
the experimental results by Wei et aZ.Q 

In the degenerate region for relatively high doping con- 
centrations, the semiclassical Mott's relationship gener- 
ally holds. P, [5^] In these cases, for a certain carrier den- 
sity (or, gate voltage), the low temperature thermopower 
scales linearly with temperature. Experimentally, the 
Mott's relationship holds very well for graphene at high 
carrier densities but breaks down close to the charge neu- 
trality point. Fig. 2 shows our results for three impu- 
rity concentrations. The linear temperature dependence 
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FIG. 1: Thermopower of clean monolayer graphene at five 
temperatures as a function of (a) carrier density and (b) chem- 
ical potential. Thermopower of impure monolayer graphene 
at four temperatures for two different impurity concentrations 
as (c) fO 12 cm" 2 and (d) 5xl0 12 cm" 2 . 
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holds even for doping close to (but still larger than) the 
peak position. For carrier densities smaller than at the 
peak position (e.g., a;~±1.4 x 1(T 4 for T=300 K and 
rii=10 12 cm -2 ) where electron-like and hole-like carri- 
ers coexist, the linear Mott's relationship is not valid. 
As impurity concentration increases, the peak position 
shifts to higher carrier densities, region of carrier den- 
sity for the violation of linear temperature dependence is 
correspondingly enlarged. 

We would like to compare our results with previous 
calculations taking into account of scattering by charged 
impurities. Hwang et al. [lj study the problem in terms 
of semiclassical Boltzmann equation approach. There, 
the high and low carrier density regions are treated sep- 
arately. Hence it is difficult to provide the dependence of 
peak position with temperature. Another work by Yan 
et aZ. [7(] starts from the microscopic Kubo's formula but 
has only focused on the low temperature limit and is also 
unable to reproduce the temperature dependence of the 
peak position observed experimentally [2tJ|- Compared 
to the above two works, we have shown that in terms 
of fully microscopic calculations incorporating the effect 
of scattering by charged impurities, both the two fea- 
tures mentioned in the beginning of this section could 
be reproduced successfully. Our results also confirm the 
picture that the peak position shifts to higher carrier den- 
sity with the increase of both temperature and impurity 
concentration. M, fijl • 



FIG. 2: Thermopower of monolayer graphene as a function of 
temperature is plotted for five electron densities and for (a) 
the clean system, (b) rii=10 12 cm -2 , (c) rii=5xl0 12 cm -2 . 



B. Unbiased bilayer graphene: The influence of 
impurity scattering 

Having achieved a reasonable success by using the 
SCBA to treat scatterings by charged impurities in a 
MLG, we shall proceed to study the bilayer system. 
We have previously studied the effect of charged im- 
purity scattering on the thermopower of gapped (bi- 
ased) BLG[ll|. In a gapped BLG, localization effect 
becomes important in the low carrier density region, 
which is out of the reach of SCBA. [54l - l56j There we re- 
stricted our attention to the systems with dilute impurity 
concentrations. 



llj However, in an unbiased BLG with a 



finite carrier density at the Fermi surface we expect to 
have a much stronger screening which allows us to con- 
sider much larger impurity concentrations. 

As a comparison to the results for MLG, we present 
in Fig. 3 thermopower of clean BLG and impure BLG 
for two typical impurity concentrations: 10 12 cm" 2 and 
5xl0 12 cm' 2 , obtained in terms of the nearest neighbor 



hopping model with 7o=3 eV and 7i=0.3 eV.[ll| For each 
impurity concentration, the carrier density dependence of 
thermopower for two typical temperatures, 300 K, 200 K 
are considered. Lower temperature requires much larger 
number of wave vectors to converge, which are difficult 
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FIG. 3: Thermopower of unbiased bilayer graphene at two 
temperatures for (a) the clean system, and the impure system 
at two different impurity concentrations as (b) 10 12 cm -2 and 
(c) 5xf0 12 cm -2 . 
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FIG. 4: Thomas-Fermi screening wave vector for an unbiased 
BLG and a MLG at 300 K and 50 K. 

observed in MLG with a smaller impurity concentration. 

The stronger screening in BLG as compared to that 
in MLG has been used to make the conjecture that 
short-range scatterers may also play an important role 
in the transport of gapless BLG in addition to the 
charged impurities. [571] It would be interesting to see how 
short-range scatterers would change the results presented 
above. 



C. Unbiased multilayer graphene 



to approach. In contrast to the MLG, impurity concen- 
tration up to 10 12 cm~ 2 does not change the qualitative 
aspects of the results. Only when the impurity concentra- 
tion is increased to be as large as 5x 10 12 cm -2 , the result 
becomes similar to that observed in experiments. fl2l fl~3|] 
In both experiments, the samples show semiconducting 
like transport behaviors with a low mobility [l2|, [III, indi- 
cating appreciable impurity concentrations and are thus 
in agreement with our results. 

The result of Fig. 3(c) is very similar to that of the 
Fig. 1(c) for a MLG except with an impurity concen- 
tration five times larger. This is due to a much better 
screening in a bilayer system. We present in Fig. 4 the 
Thomas-Fermi screening wave vectors of both systems at 
two different temperatures: 300 K and 50 K. At 300 K, 
qrF of BLG is more than five times that of MLG around 
the charge neutral point. This difference becomes even 
larger at 50 K. For the carrier density where the ther- 
mopower peaks (x ^±2xl0~ 4 ) for BLG, qxF of BLG 
is still more than two times that of MLG. According to 
Eq. (16), Eq. (17), and Eq. (18), the effective strength 
of impurity scattering is approximately proportional to 
niq^p. Hence a BLG with a much larger concentration of 
charged impurity shows the same qualitative behavior as 



Now, we study the evolution of thermopower as a func- 
tion of layer number for an unbiased multilayer graphene. 
From now on we would only consider clean graphene mul- 
tilayers. First, we consider the simple nearest-neighbor- 
hopping model, ignoring temporarily the screening ef- 
fect. As shown in previous works, this simplified multi- 
layer graphene model could be decomposed into indepen- 
dent subsystems of MLG and BLG.[30|, [H M, Szl The 
two linear response coefficients L\\ and L 12 are then ob- 
tained as summations of and which are the 
corres pon ding values for the various subsystems labeled 
by m.[30[ Thermopower of multilayers is then obtained 
from Eq. ([9]). The room temperature thermopowers ob- 
tained in this way are shown in Fig. 5(a) with layer num- 
ber up to six. Peak value of the thermopower does not 
show very large variation with layer number. Different 
from MLG, thermopower for larger layer number sam- 
ples usually show secondary peaks associated with the 
onset of higher conduction or valence bands contributing 
to transport. fill] 

Results presented in Fig. 5(a) does not take the ef- 
fect of screening into account, so electrons on different 
layers feel the same electrostatic potential. However, as 
discussed in Sec. II C, multilayer graphene with layer 
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80 




FIG. 5: Room temperature (300 K) thermopower of clean 
unbiased multilayer graphene. (a)Nearest neighbor hopping 
model without including the screening effect or charge redis- 
tribution between layers. (b)Nearest neighbor hopping model 
with the screening effect included. (c)Full SWMcC model 
with screening effect incorporated. 



number larger than two should have different carrier den- 
sities for different layers. When this effect is treated 
self-consistently, results for the multilayer graphene de- 
scribed in terms of the simplified model are presented in 
Fig. 5(b). Except for some quantitative difference, the 
qualitative behavior is the same as in Fig. 5(a). The 
positions and magnitudes of various peaks are separated 
into two groups depending on whether N is even or odd. 
For both cases, the peak value of thermopower increases 
and the peak position continuously shifts to higher carrier 
density with increasing layer number. In the presence of 
screening effect, site energies of the various layers become 
unequal, hence the subsystem decomposition is not valid 
and Fig. 5(b) is obtained starting from the full model. 

Fig. 5(c) shows the thermopower of multilayer 
graphene described in terms of the more general 
Slonczewski-Weiss-McClure model [UM H, Em [Hj], 



with the screening effect incorporated. Since the particle- 
hole symmetry is now explicitly broken, the curve be- 
comes asymmetric with respect to the point of zero car- 
rier density. Specifically, the positive and negative peak 
value becomes unequal and the position of zero ther- 
mopower is shifted away from zero doping. Another 
qualitative difference as compared to the results obtained 
from the simplified model is that, the peak value of the 
even-layer (N=2, 4, 6) and odd-layer (iV=3, 5) systems 
show opposite behaviors compared to Fig. 5(b). This 
nontrivial result arising from further neighbor hoppings 
contained in the SWMcC model needs experimental con- 
firmation. 

The reason for the above grouping between even and 
odd layer systems could be understood qualitatively from 
the subsystem decomposition in Sec. IIB. Thermopower 
of various bilayer like graphene subsystems form a se- 
ries which show monotonous dependency on the effec- 
tive interlayer hopping 71 m in Eq. ([5]). Thermopower of 
monolayer like subsystem does not fall into this series. 
Since monolayer like subsystems are contained only in 
odd-layer graphene, it makes sense that the even layer 
and odd layer multilayer graphene should group sepa- 
rately. Though self-consistent calculation and inclusion 
of further neighbor hopping both make the subsystem de- 
composition invalid, numerical results contained in Figs. 
5(b) and 5(c) are still consistent with the above picture. 



D. Biased multilayer graphene: The electronic 
structure 

Transport properties, like conductivity and ther- 
mopower, sensitively depend on the underlying electronic 
structure. In order to understand the thermopower of 
biased multilayer graphene to be presented in the next 
section, we first discuss in this section the electronic 
structure of a biased multilayer graphene, taking tri- 
layer and quad-layer graphene systems as two exam- 
ples. In the presence of an electric field perpendicular 
to the layer plane, inversion symmetry of a multilayer 
graphene is broken explicitly. In this case, even for the 
bilayer system, charge redistributes to screen the elec- 
tric field. As in the above section, we take this screening 
into account in terms of the Hartree approximation by 
treating separate layers as parallel plates with certain 
density of net charge which are to be determined self- 
consistently. [H-il m HI 

We focus on two related questions. One is the robust- 
ness of the Dirac-like linear dispersive band in the pres- 
ence of an electric field for a multilayer graphene with 
an odd number of layers. The other is the possibility of 
opening a full ener gy gap by an external electric field, 
similar to the BLG.|lllll [l| 

We find that the effect of temperature on the electronic 
structure is very tiny. Thus, only zero temperature re- 
sults are presented here. For the biased systems, without 
loss of generality, the bare interlayer potential energy dif- 
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FIG. 6: Low energy band structures for unbiased and bi- 
ased trilayer graphene. (a) Vb=0, n=0, /i~-20.8 meV. (b) 
V o =0, n=2xl0 13 cm" 2 , /x~542.1 meV. (c) V =l eV, n=0, 
/i~2.4 meV. (d) V =l eV, n=2xl0 12 cm -2 , /x~45.8 meV. n 
is the average carrier density per layer. TK and KM denote 
respectively segments along the high symmetry lines in the 
two dimensional BZ connecting two of three highly symmet- 
ric points F, K or M. 



ference induced by the electric field is taken as Vo=l eV. 

For an unbiased charge neutral trilayer graphene, as 
shown in Figs. 6(a), the system is gapless. But when the 
carrier density is increased to be as high as 2xl0 13 cm -2 , 
a small gap about 5 meV is induced below the chemical 
potential. Two Dirac cone like structures still exist in 
the system, but they do not touch at the cone vertex and 
are usually immersed in bilayer like bands. Applying a 
perpendicular electric field as high as Vq=1 eV, a small 
gap about 17 meV is induced even for the neutral system 
x=0, as shown in Figs. 6(c). The magnitude of the gap 
enhances slightly with the increase of x as shown in Fig. 
6(d). 

For the biased quad-layer graphene, as shown in Figs. 
7(c) and 7(d), a full gap is not observed. For the un- 
biased quad-layer graphene, though the neutral system 
is still gapless (Fig. 7(a)), a full gap below chemical 
potential is observed when the carrier density is high 
enough, as shown in Fig. 7(b). The latter feature, to- 
gether with that in Fig. 6(b), arises from the screening 
effect. [3l|, HH, IH, HI] Screening effect for the unbiased 
multilayer graphene systems for N>2 matters because, 
for layer number larger than two, N/2 (for even N) or 
(N + l)/2 (for odd N) non-equivalent layer sets appear 
which usually have different on site energies and electron 
densities. 

Finally, we present in Fig. 8 the density of states 
(DOS) of the trilayer and quad-layer graphene for several 
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TK K KM TK K KM 



FIG. 7: Low energy band structures for unbiased and biased 
quad-layer graphene. (a) Vo=0, n=0, /i~-20.8 meV. (b) Vo=0, 
n=2.5xl0 12 cm" 2 , /x~124.3 meV. (c) Vb=l eV, n=0, M ~15.5 
meV. (d) V =l eV, n=2.5xl0 12 cm" 2 , ^~43 meV. Meanings 
of FK and KM are the same as in Fig. 6. 

typical parameter sets. The corresponding room tem- 
perature chemical potentials are marked by the dotted 
vertical lines. As compared to the unbiased systems, ex- 
ternal bias introduces some Van Hove singularities in the 
density of states. 

E. Biased multilayer graphene: The thermopower 

Previously, we have shown that the opening of a gap 
in biased BLG greatly enhances thermopower of the 
svstem flll |59| . Calculations in the previous section show 
that a small gap opens also in other multilayer graphene 
systems, such as the trilayer graphene. It is thus inter- 
esting to see whether or not the thermopower is likewise 
enhanced in these systems. 

Results for multilayer graphene systems at room tem- 
perature with layer number up to 5 are presented in Fig. 
9. Strength of the external electric field is taken as Vq—1 
eV. Fig. 9(a) contains results obtained from the nearest 
neighbor hopping model. Fig. 9(b) is for the SWMcC 
model. 

Compared with zero bias results, the large enhance- 
ment of thermopower under a bias is specific to the BLG 
system. Though peak value of thermopower is slightly 
enhanced in trilayer graphene (Fig. 9(c)), it decreases 
under a bias for both the quad-layer and quintuple-layer 
systems. 

In a biased BLG, peak value of room temperature ther- 
mopower is roughly proportional to the gap size. An ex- 
ternal field as large as Vq=1 eV opens a gap of about 
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FIG. 8: Low energy density of states (DOS) for unbiased 
and biased trilayer and quad-layer graphene, normalized to 
be corresponding to 6 (8) carbon atoms for trilayer (quad- 
layer) graphene per spin. (a)Unbiased trilayer graphene, 
with n=2xl0 12 cm -2 . Room temperature chemical po- 
tential is ^t(300K)~ 101 meV. (b)Biased trilayer graphene, 
with n=2xl0 12 cm -2 and Vb=l eV. m(300K)~ 33.2 meV. 
(c)Unbiased quad-layer graphene, with n=2.5xl0 12 cm -2 . 
/i(300K)~ 121 meV. (d)Biased quad-layer graphene, with 
n=2.5xl0 12 cm" 2 and V =l eV. ^(300K)~ 42.3 meV. Dotted 
vertical lines mark position of the corresponding room tem- 
perature chemical potentials. Insets are enlargements of the 
small energy regions. 



288 meV, which yields a peak room temperature ther- 
mopower of about 412 /uV/K.|Tl} From the previous sub- 
section, Vb=l eV gives a gap of about 20 meV for trilayer 
graphene. Following this thread of argument, room tem- 
perature thermopower of the biased (Vb=l eV) trilayer 
graphene would be around 50 /iV/K, which is in qualita- 
tive agreement with Fig. 9(b). 

For layer number larger than three (e.g., Fig. 7(c) 
and 7(d) for the quad-layer system), usually no full gap 
opens in the biased multilayer graphene system since the 
bottom and top layers almost decouple for so thick sam- 
ples. In addition, as could be seen from Figs. (6), (7) 
and (8), the chemical potential is usually not situated in- 
side the gap, even though a gap opens for certain doping 
and external electric field. Both of these effects suppress 
the thermopower for biased multilayer graphene systems 
with iV>3. 

To confirm the self-consistency of the results, we an- 
alyze the applicability of the Mott's formula, which re- 
lates the thermopower to the longitudinal conductivity 
as 0, SHU 
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FIG. 9: Room temperature (300 K) thermopower of bi- 
ased multilayer graphene as a function of carrier density, 
(a) Nearest neighbor hopping model with the screening effect 
taken into account. (b)SWMcC model with screening effect 
incorporated. (c)A trilayer graphene with and without a bias 
potential. 



The conductivity is obtained from the linear response 
coefficient as a—e 2 Ln/T. We present in Fig. 10 re- 
sults obtained by a full microscopic calculation and by 
fitting the Mott's formula, taking the room temperature 
thermopower of biased (Vq=1 eV) trilayer and quad-layer 
systems as two examples. The Mott's formula holds well 
at temperatures smaller compared to the Fermi energy. 
As shown in Fig. 10, the qualitative behavior of ther- 
mopower is reproduced very well by the Mott's formula 
even at the room temperature. Only close to the peak 
position, is the difference appreciable. 

From experiences gained on the monolayer and bilayer 
graphene, it is safe to say that impurity scattering would 
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FIG. 10: Room temperature (300 K) thermopower of biased 
(Vq — 1 eV) (a) trilayer and (b) quad-layer graphene. A com- 
parison of results from full microscopic calculation (labeled 
as "Kubo") and a fitting from the Mott's formula (labeled as 
"Mott") is made. The SWMcC model including the screening 
effect is considered here. 



not enhance thermopower much. Thus, we conclude that 
biased BLG systems has the largest room temperature 
thermopower among all the multilayer graphene systems. 



IV. SUMMARY 

In this work, we systematically calculate thermopower 
of biased and unbiased multilayer grphene systems. The 
effect of screening of an external electric field is taken 
into account self-consistcntly under the Hartree approxi- 
mation so that charge densities are different between in- 
ner and outer layers. Both the model with only near- 
est neighbor hopping and the more general SWMcC 
model with further neighbor hoppings are considered. 
The effect of impurity scattering is considered for mono- 
layer and unbiased bilayer graphene in terms of the self- 
consistent Born approximation. For monolayer graphene, 
only when the effect of impurity scattering is taken into 
account, we could obtain results consistent with experi- 
ments. The electronic structure and thermopower of bi- 
ased multilayer graphene systems are calculated in which 
the screening effect is self-consistently incorporated. The 
biased bilayer graphene shows the largest room temper- 
ature thermopower. 
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